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We perform full-MHD simulations on various initially helical configurations and show that they 
reconfigure into a state where the magnetic field lines span nested toroidal surfaces. This relaxed 
configuration is not a Taylor state, as is often assumed for relaxing plasma, but a state where the 
Lorentz force is balanced by the hydrostatic pressure, which is lowest on the central ring of the 
nested tori. Furthermore, the structure is characterized by a spatially slowly varying rotational 
transform, which leads to the formation of a few magnetic islands at rational surfaces. We then 
obtain analytic expressions that approximate the global structure of the quasi-stable linked and 
knotted plasma configurations that emerge, using maps from S^ to S‘^ of which the Hopf fibration is 
a special case. The knotted plasma configurations have a highly localized magnetic energy density 
and retain their structure on time scales much longer than the Alfvenic time scale. 


Understanding the types of structures in magnetic 
fields that occur in magnetohydrodynamics (MHD) is of 
fundamental importance for nuclear fusion mia and as¬ 
trophysics (US]. Helicity-constrained, unbounded exci¬ 
tations in plasmas are present in a wide range of scales, 
from underdense bubbles emitted from active galactic nu¬ 
clei ( r\j 10 kpc) [7], through magnetic structures ejected 
from the solar corona 10^“^ km) [8] to the struc¬ 
ture in fusion reactors such as the spheromak [9] and 
field-reversed configurations nni m), and the plas- 
moids in dense plasma focus (DBF) experiments mm) 
HD. There exist many analytical solutions for the field 
in toroidal confinement vessels [12] and bounded domains 
m, and even confinement vessels in the shape of a knot 
El- There are also analytical expressions for unbounded 
force-free fields [T5|, but no analytical expression has 
been found for a localized field that agrees with observed 
structures seen in unbounded plasmas. 

Magnetic helicity, defined as = / A • B d^x, where 
A and B are the vector potential respectively the mag¬ 
netic field, was recognized by Woltjer to be an invariant 
of an ideal plasma [16]. The identification of helicity 
as linking of magnetic field lines by Moffatt mi gave a 
clear topological interpretation. Given the topological 
nature of this invariant, Kamchatnov used the structure 
of the Hopf fibration to construct a topological soliton 
in ideal MHD [18]. Recently this work was generalized 
by Thompson et al HSj. This structure has not been 
described in resistive MHD, but also there helicity and 
magnetic topology play an important role in constraining 
magnetic relaxation naiinHisi. In order to understand 
the effect of helicity in resistive plasmas we simulate the 
time evolution of various helical initial conditions and 
find that each of them evolves towards an ordered state 
of nested toroidal magnetic surfaces. 

We simulate the plasma dynamics using the PENCIL 
CODE [24]. With this code we solve the resistive MHD 
equations in dimensionless form for an isothermal plasma 


in a fully periodic box of volume (27r/o)^ (see supple¬ 
ment). We choose as initial conditions simple configura¬ 
tions that are clear examples of fields containing helicity; 
rings of flux that are all linked and/or twisted. We start 
simulations with n identical magnetic flux tubes that are 
all linked, with n ranging from 2-6. The flux tubes have 
magnetic field of IBq at the center of the tube, a radius 
of V 2 I 0 , and a Gaussian intensity profile with charac¬ 
teristic width of 0.16/o- For the n = 3 configuration 
we also vary the twist T which indicates the number of 
windings of a field line around the center of the tube 
as it passes around the tube once, further increasing he¬ 
licity. Two initial conditions are shown in figure (a) 
and (c). The velocity is initially zero everywhere and 
the density p is set to 1 uniformly in the initial condi¬ 
tion. The kinematic viscosity and magnetic diffusivity 
are 2 x 10“^, giving a magnetic Prandtl number of unity. 
The magnetic helicity of the initial conditions is given by 
i7m = — n)4>^ + nT4>^, where 4> is the magnetic flux 

through a single ring (see supplement). 

The configurations evolve in a similar fashion, which 
can be divided into two regimes, reconnection and resis¬ 
tive decay, as shown in figure [^e)-(f). We use Alfvenic 
time tA = l/(2\/27r/o), scaled by the length of a flux tube. 
The tubes first contract, as this lowers the magnetic en¬ 
ergy. This process is further detailed in the supplement. 
The higher the initial helicity, the less energy can be lost 
through reconfiguration. Figure (g) shows the evolu¬ 
tion of several related quantities for the simulation with 
n = 3 and T = 1.8. 

In order to analyze the emerging plasma configuration 
we take a detailed look at the simulation with n = 3 
and T = 1.8 at time t = 22.5tA- The magnetic energy 
is highly localized (figure]^ (a)), falling off from the cen¬ 
ter. Remarkably, from the chaotic collapse of the initial 
condition, containing only a discrete rotational symmetry 
around the z-axis, an ordered magnetic structure emerges 
that is roughly axisymmetric and where field lines span 
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FIG. 1. Simulated configurations and evolution of magnetic fields. (a),(b) Initial condition and state at t = 22.5tA for the 
n = 3 and T = 1.8 simulation, (c)-(d) The same for the n = 6 and T = 0 simulation. In (a) and (c) a magnetic isosurface of 
|B| = O.li^o is shown to indicate the boundary of the flux tube. In (b), (d) the lines represent the magnetic field lines. The 
outer field lines are partially transparent to not obstruct view of the central configuration, (e) Decay of magnetic energy for 
the simulations with T = 0 and n ranging from 1-6, and (f) the simulations with n = 3 and T ranging from 0-4.4. The shaded 
region indicates reconfiguration, after that resistive decay takes over, (g) The evolution of the average of magnetic energy 
(B^)/(Bo), normalized helicity (A-B)/(B^), helicity (A'B)/(Bo) and velocity (|v|) for the simulation with n = 3 and T = 1.8. 


invariant tori. These are toroidal surfaces spanned by 
magnetic field lines and are often described in the con¬ 
text of toroidal fusion devices [25] . Four toroidal surfaces 
are shown in figure(a). 

With higher initial helicity this structure appears 
sooner and is more pronounced. Invariant tori are ob¬ 
served in all simulations except the n = 2 simulation 
which was stopped at t = GO^a- In the n = 3 and 
T = 0 simulation tori were found only after t = 54tA, 
but in all other simulations this structure appears before 
t = 22.5tA and remains. Invariant tori are also observed 
in simulations using different helical initial conditions, 
such as a single twisted ring and a trefoil knotted flux 
tube (see supplement). 

The initial reconfiguration of the rings induces pressure 
waves traveling through the periodic simulation volume. 
However, these pressure waves do not significantly affect 
the magnetic structure. To investigate the role of pres¬ 
sure in the simulation we average out these waves by aver¬ 
aging 365 snapshots between t = 27.5tA and t = 3b.St a- 


Figure |^(c) shows the averaged pressure, which is lowest 
on the magnetic axis of the structure. An ambient pres¬ 
sure Poo is therefore inherent to the structure. The force 
due to the pressure gradient is balanced by the Lorentz 
force, which makes the structure quasi stable. In figure 
1^ (a) and (b) we show the average radial component of 
the Lorentz F£, and minus the average radial component 
of the pressure gradient — VP^ ia the x,p-plane pass¬ 
ing through the center of the structure (top view of the 
torus). 

Note that the lowered pressure in the structure is con¬ 
sistent with the virial theorem [221 EH that states that a 
free plasma cannot uphold an inerease in pressure solely 
by internal hydromagnetic forces. The region of highest 
magnetic field strength is near the geometrical center of 
the tori, where the pressure is unchanged. 

The balance of magnetic and hydrostatic forces indi¬ 
cates that the magnetic field forms self-stable, localized 
structures in MHD equilibrium with ambient pressure 
Poo- These equilibria feature rich dynamics such as the 
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FIG. 2. The simulation with n = 3 and T = 1.8 at time 
t = 22.5tA. (a) The magnetic field contains invariant tori. 
Every surface is a single integral curve of the field of length 
lOOO/o colored differently for clarity. The surfaces are clipped 
to show the nested configuration, (b) Surfaces of constant 
magnetic field strength. A single torus is shown in black to 
indicate the scale of the magnetic structure. 



FIG. 3. (a) The radial component of minus the gradient of 
the averaged pressure field, and (b) the radial component of 
the averaged Lorentz force, taken in the x, y-plane (top view). 
The geometrical center of the tori is taken as the origin , and 
is marked by the blue dot. (c) The pressure field in the x^z 
half-plane, showing a lowered pressure in the center of the 
magnetic surfaces. 


formation of magnetic islands at rational surfaces that 
are an area of intense research [25l [28] . 

In order to investigate the nature of this equilibrium we 
construct a Poincare plot of the field in figure As seed 
points we choose 500 points on a line from the geometrical 
center of the tori and through the magnetic axis, starting 
on the magnetic axis and moving outward. We label this 
direction x*, and the direction perpendicular to that and 
out of the plane of the torus we call z*. The field lines are 
traced for a distance of 4000/o, and the positions where 
they cross the plane defined by x* and 2 ;* are marked. 
The Poincare plot is shown in figure(a). We show the 
rotational transform i [29] of the corresponding field line 
in figure 1^ (b) (see supplement for calculation). 


As expected from [25], where the rotational trans¬ 
form crosses rational values we observe magnetic islands. 
Lines are drawn indicating where the rotational trans¬ 
form crosses the values 8/9, 7/8, 6/7 and 5/6. As ex¬ 
pected the number of islands observed is equal to the 
denominator of the rotational transform. Even though i 
crosses a few rational surfaces, the value varies less than 
10%. In a tokamak equilibrium, where the inverse of 2 , 
the safety factor is used, this value typically varies 
much more [30] . We note that the fact that our pressure 
plots result from an averaging over time implies that we 
cannot resolve the fine structure in the pressure, such as 
possible discontinuities in pressure over specific irrational 
KAM surfaces as described in m- 

The magnetic field strength and the x*, ^*, and 2 :*- 
components (with perpendicular to 2 :* and x*) at the 
position of the seed points are shown in figurej^ The field 
varies continuously over the surfaces, and the magnitude 
is higest in the geometrical center of the structure. 

The part of the magnetic field that forms toroidal mag¬ 
netic surfaces is reasonably axisymmetric, and could in 
principle be approached by a solution to the Grad Shafra- 
nov equation EH- This would however not capture the 
large part of the field outside of this ordered region. In¬ 
stead we want to point out a curious resemblance be¬ 
tween the structure of the Hopf fibrations and the fields 
observed here. 

Non-null-homotopic maps (functions) from (hyper¬ 
sphere) to 5'^ (sphere) such as the Hopf map [32] feature 
a topolocical structure resembling the observed plasma 
sctructures. The fibers (pre-images of points on S^) of 
the map are continuous curves that lie on the surfaces of 
nested tori. Furthermore, every fiber is linked with every 
other one, with linking number depending on the map. 
Through stereographic projection from to the fiber 
structure of this map can be translated to a vector field 
in whose integral curves are the fibers of this map, 
(derivation in the supplement). Moreover, the obtained 
field is smooth, continuous, divergenceless, has helicity, 
and the field lines lie on the surfaces of nested tori. 

This curious structure was used by Kamchatnov to de¬ 
scribe a soliton in ideal MHD, where the fluid velocity 
is parallel to the magnetic field everywhere m- Inde¬ 
pendently, Rahada, used the structure of the Hopf map 
to construct full radiative solutions of Maxwell’s equa¬ 
tions [33] [34]. Kamchatnov’s solution was generalized by 
Sagdeev [35], and a similar extension of Rahada’s fields 
was described by Arrayas and Trueba [36] . 

The analytical form of this vector field in is given 
by: 


4r^Va 

7r(rg +r2)3 


( 2{L02roy - ujixz) \ 
-2{w2rox + ojiyz) 

M-rl+x^ + y^-z^)J 


( 1 ) 


This field is cylindrically symmetric around the z-axis 
(see supplement). It has a finite magnetic energy, as can 
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FIG. 4. Poincare plot and properties of the force-balanced 
toroidal structure, (a) Poincare plot of the magnetic held, 
(b) The value of the rotational transform for each magnetic 
surface. Rational values are indicated by the labeled hori¬ 
zontal lines, and the positions where they cross are shown 
by vertical lines, (c) Value of the magnetic held strength, 
and the components at each position, (d) The magnetic held 
strength and components of the magnetic held for the analyt¬ 
ical expression of a held with the same energy and rotational 
transform. 


bee seen from: 


J d^x = arQ{ujl + ujI), ( 2 ) 

and nonzero helicity, given by: 

-ffm = arQUiL02, (3) 

and like the held in the simulation it tends to zero away 
from the center. 


In our observed structure, the huid velocity is neither 
parallel nor proportional to the magnetic held, making 
this structure fundamentally different from the structure 
Kamchatnov described [18]. Nevertheless, the magnetic 
topology of held lines lying on nested toroidal surfaces, 
the magnetic energy localized in the center, the near con¬ 
stant rotational transform, and the direction of the mag¬ 
netic held, even outside of the area that forms magnetic 
surfaces, all agree qualitatively with the toroidal struc¬ 
ture described by equation To quantify this claim 
we extract from the simulation the parameters cji, UJ 2 , 
and ro, needed for equation 16 (method described in the 
supplement), and show that there is overall agreement. 

For the simulation with n = 3 and T = 1.8 this yields 
values of ro = 0.78, uJi = 0.24 and uj 2 = 0.27. Parameters 
for the other simulations are quite similar (see supple¬ 
ment). The analytical magnetic held is shown in hgure 
1^ (d) for the same positions as the extracted held in (c). 
Even though there are differences in the magnitude of the 
components, there is broad agreement, which is quite re¬ 
markable for a routine that only uses three independent 
variables that are not ht, but calculated from select pa¬ 
rameters extracted from the simulation. As time elapses 
ro increases and ujiIuj 2 decreases. This change is such 
that over 45tA ^o increases by 35% and co’i/co ’2 decreases 
by 50%. 


We have shown that reconnection of helical helds in re¬ 
sistive MHD causes the emergence of a self-stable toroidal 
magnetic held in force equilibrium. This equilibrium re¬ 
sults from a balancing of magnetic forces and the pres¬ 
sure gradient, and has a minimum in pressure on the 
magnetic axis. Note that this is not a Taylor state, and 
the pressure prohle is inverse to the pressure enforced in 
a Tokamak reactor. In the quasistable state there is rich 
dynamics such as the emergence of magnetic islands at 
rational surfaces. 

Furthermore, we obtained an analytic expression for a 
magnetic held whose held lines lie on nested tori, requir¬ 
ing only three independent parameters. This held is a 
good approximation for the plasma conhgurations that 
emerge in the simulations, where a signihcant portion 
of the magnetic held lines reconhgure to lie on nested 
toroidal magnetic surfaces. We have observed the forma¬ 
tion of this self-stable structure for various initial plasma 
conhgurations containing helicity. This indicates that 
this structure is a fundamental self-conhning conhgura- 
tion that we predict to occur in situations where there is 
unbounded plasma containing helicity. 
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APPENDIX: SUPPLEMENTAL MATERIAL 
MAGNETIC HELICITY 

We calculate the helicity of the initial condition con¬ 
sisting of n linked rings each with a twist T, each carrying 
the same magnetic flux This calculation holds for thin 
rings. The helicity integral 


plasma in sixth order finite differences in space and third 
order in time. We solve the equations in a fully peri¬ 
odic box of size (27r/o)^ with 256^ meshpoints. We solve 
for the vector potential A instead of the magnetic field 
B, and use the resistive gauge. For an isothermal gas 
the pressure is given hy p = pc^ , where Cg is the sound 
speed, set to unity and p is the density. The equations 
we solve are: 


[ A • B d^x (4) 

Jv 

can be split into the contributions of every single ring 
separately, since the magnetic field is zero outside the flux 
rings. For a thin ring the integral can be approximated 
by 

^ A • dl. (5) 

Here 4> is the magnetic flux in a single ring and the in¬ 
tegral is evaluated along a closed curve running along 
the center of the flux ring. Through Stokes’ theorem the 
integral can be rewritten to: 

^ V X A • da = $ ($encl) , (6) 

where the definition for the vector potential V x A = 
B has been used, and 4>enci indicates the magnetic flux 
tbrought the surface enclosed by the curve. This has two 
contributions 


—-V X B — A = 0, 

ot 


(9) 


Dv 


+ c^Vlnp - 


-(j X B) - -Fvisc = 0, 
P P 


( 10 ) 


D Inp 
Dt 


+ V-v = 0. 


( 11 ) 


Here v is the plasma fluid velocity, t is time and the 
operator ^ + v • V) indicates the material deriva¬ 

tive. The viscous forces are given by Fvisc = V • 2i'pS 
with v the kinematic viscosity and S the traceless rate of 
strain tensor (S = l/2{vi^j-\-Vj^i) — l/‘^5ijV'w). The equa¬ 
tions are solved on a grid size 256 ^. The gridsize of 256 ^ 
is more than adequate to resolve the spatial structures 
that we are interested in. 


OPEN BOUNDARY CONDITIONS 


^encl — ^self ^others5 (^) 

where 4>seif indicates the contribution of flux by the field 
of the ring itself, caused by twist in the ring, and 4>others 
indicates flux caused by other flux tubes passing through 
the center of the ring integrated over. 

The contribution of 4>others is proportional to the num¬ 
ber of rings passing through the ring in question, and is 
equal to 4>others = (n —1)4>. The contribution 4>seif of self¬ 
linking in a twisted flux tube is proportional to the twist; 
^seif = F4>. This defines the helicity of a single ring, the 
helicity of the entire configuration is simply found by 
multiplying by n: 

= (n^ - n)$2 + (8) 

The parameters n and T increase the amount of mag- 
netic helicity in the initial condition. For thin rings T 
only changes the magnetic helicity of the initial configu¬ 
ration, not the magnetic energy. 

EQUATIONS SOLVED BY THE PENCIL-CODE 

The PENCIL-CODE is used to solve the resistive mag- 
netohydrodynamical (MHD) equations for an isothermal 


To assess the effect of the boundaries on the observed 
configuration we perform simulations using open bound¬ 
ary conditions in addition to the simulations with peri¬ 
odic boundary conditions reported in the main article. 
We simulated the n = 3, T = 1.8 initial condition using 
open boundary conditions. These conditions are imple¬ 
mented by imposing the condition that the fields (velocity 
V and magnetic field B) are perpendicular to the bound¬ 
ary, allowing magnetic flux to exit the simulation volume, 
and a smooth variation of the density across the bound¬ 
ary by setting the first derivative across the boundaries 
to zero. 

Eigurej^ shows the time evolution of relevant parame¬ 
ters. There is almost no difference between open bound¬ 
aries and the periodic boundaries used, with the notable 
exception of the velocity field. The velocity field in the 
periodic boundaries is caused by the fluid motion in the 
pressure waves that are created in the initial collapse of 
the rings, but in the open boundary simulations these 
waves escape the volume. 

Also the magnetic fields are very similar in topology. 
The inset in figureshows the magnetic field and toroidal 
surfaces, using the same seed points for the stream trac¬ 
ing as in figure 2 (a) of the paper. 
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FIG. 5. Comparison of periodic boundary conditions to open 
boundary conditions, (a) Quantities in the simulations decay 
in a similar fashion. The cyan fat curves are from the peri¬ 
odic boundary conditions simulation, whereas the thin black 
curves are from the open boundary conditions simulation, (b) 
Field structure in the open boundary conditions simulation at 
t = 42tA- The same nested structure appears. 


TREFOIL KNOT AND SINGLE TWISTED RING 

In order to assess the generality of the observed emerg¬ 
ing plasma configuration we also performed two other 
simulations with initial conditions that have helicity: The 
trefoil knot and a single twisted torus. The trefoil is 
parametrized as described in reference [22] in the main 
paper. The field is set to constant ISq inside the rings 
and zero outside. The twisted torus is parametrized in 
the same way as the twisted rings, but the ring is set in 
the x,^-plane. The twist is set to T = 3.1. The initial 
conditions can be seen in figure]^ (a) and (c). All the 
simulation parameters are identical to the ones used in 
the paper: fully periodic simulation volume of (27r/o), 
isothermal gas with viscosity and magnetic diffusivity 
of77 = z^ = 2xlO“^. For easier comparison the same 
alfvenic timescale is used. 

The fields evolve in the same manner as the linked 
fields rings in the paper, with an initial fast drop in mag¬ 
netic energy, followed by a much slower decay. They also 
show the self-formed nested toroidal surfaces that ap¬ 
peared in the ring simulations, as is shown in figure]^ (b) 
for the trefoil simulation (at t = 45tA) and (d) for the 
twisted torus simulation (at t = 22.5tA)- 


VIDEO 

The supplemental material contain a video visualizing 
the time evolution of the simulation with n = 3 and 
T = 3.5. This video shows the collapse of the rings and 
the emergence of the stable structure. Table [I| shows at 
what times which processes are shown in the movie. 



FIG. 6. Initial conditions, final states, and evolution of mag¬ 
netic energy for the trefoil and twisted torus simulation, (a). 
Initial state of the trefoil simulation. The field strength is 1 
inside the tubes, and zero outside, (b). Final state of the tre¬ 
foil simulation, showing the nested toroidal surfaces spanned 
by the field lines (clipped to show nesting) and 70 other field 
lines colored by magnetic intensity, (c) Initial condition for 
the twisted torus simulation. The field has a Gaussian pro¬ 
file. An isosurface of magnetic field strength 0.1 shows the 
approximate edge of the flux tube. The field strength is 1 in 
the center of the flux rings. A line of 100 closely spaced field 
lines shows the twist of the field, (d). Final configuration of 
the twisted torus simulation, visualized in the same way as 
(b). (e). The decay of magnetic energy for the n — 3^T — 1 
simulation, the trefoil simulation, and the twisted torus sim¬ 
ulation. 


TIME EVOLUTION OF FIELDS 

Figure shows in more detail the time evolution of 
different initial conditions, including the evolution of just 
a single ring. The first row shows the initial condition. In 
















FIG. 7. Time evolution of the magnetic fields. The properties of the initial condition are given on the top of each column. 
Time increases in the downward direction. In the initial condition an isosurface for the magnetic energy is drawn indicating 
the extent of the ring. The configurations with more than one ring show the emergence of a localized linked field configuration, 
the single ring does not. At t = 22.5tA details of the initial structure are lost and a universal structure of linked field lines on 
nested toroidal surfaces appears. 
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TABLE L Video events and time, made from the simulation with n = 3 and T = 3.5 


time (movie) 

time (simulation) 

event 

0-5s 

0-4.5tA 

initial collapse (slowed down) 

5-14s 

4.5-34.7tA 

formation of structure 

14-17S 

34.7tA 

plane showing B 

17-20S 

34.7tA 

zoom into center of structure 

20-26S 

34.7tA 

held line spanning hrst toroidal surface 

26-32S 

34.7tA 

held line spanning second toroidal surface 

32-37S 

34.7tA 

held line spanning third toroidal surface 


the second row, taken at t = 0.5tA it is seen that the rings 
have contracted. This is to be expected, as this lowers 
the magnetic energy. Another way of understanding this 
is by looking at the currents that must source this field. 
A flux tube is sourced by ring-currents that run around 
the minor radius of the tube. These currents run parallel, 
and therefore must attract one another. The net effect is 
a force that causes the magnetic flux tube to contract. 

If one looks closely at the difference between the n = 1 
and n = 2 situation one can see the physical manifesta¬ 
tion of helicity conservation in action. The single ring 
can contract onto itself, but the two rings cannot cross 
each other. 

When the rings have collapsed, shown in the third row, 
the single ring can then lose a large amount of field, that 
is carried away from the center by fluid motion. This 
cannot happen in the linked rings. 

The fourth row shows the state at time t = 4.5tA- It 
is clear that in the linked configuration the energy is still 
confined in the center, but the initial structure is still 
seen. At time 22.5tA most of the initial structure is lost, 
and only a localized, helical configuration is seen for the 
linked configurations, whereas little structure is visible in 
the simulation with a single ring. 


CALCULATION OF THE MAGNETIC FIELD 


Consider an adaptation of the Hopf map 


._ ^(2) ^ 



( 12 ) 


Here indicates the operation on a complex number 
such that only its argument is multiplied by the real num¬ 
ber w. The two complex numbers ( 2 ) 1 , 2 ) 2 ) C with 
kil + kil = 1 define a point on the three-sphere 
and : C ^ >5^ is inverse stereographic projection 

from the north pole. From this definition follows that 
(^zi^ Z 2 ) = for all values 

of 0. Consequently the pre-image of a point on S‘^ is a 
continuous curve in . Furthermore, if uji = UJ 2 each 
curve is a great circle in S'^, and the map is in essence 
identical to the Hopf map. For uji ^ uj 2 the pre-image of 
a point in S‘^ is a curve in that oscillates in the zi- 


direction with frequency uoi and in the 2 ) 2 -direction with 
frequency Co’ 2 . 

To construct a field in define the following map: 

0 : ^ C, (/) = o o ^(3)-^^ ^^ 3 ) 


with ^ ^ 5'^ inverse stereographic projection 

(also from the north pole of S^) to the three-sphere. This 
map has the following explicit form: 


(l){x,y,z) 


{2{x + 

(2z -h i(r^ — 


(14) 


where = x^ -j- y^ -h z^. We get a vector field in by 
calculating 


^ V(/) X Vcj)"" 

27ri (1 -h 00*)^ 


(15) 


Here ^/a is a constant such that the magnetic field has 
correct dimensions. The cross product between Vcj) and 
V0* is a vector in the direction that (j) remains constant. 
The curves of constant 0 are continuous, oscillating and 
closed (or, for incommensurate uJi and Co’ 2 , dense in a 
compact subspace of S^) curves in 5'^, and will thus be 
so too in U 00 . We calculate this field, and then scale 
it with a factor tq by substituting {x,y^z) ^ 7 ^) 

to obtain the expression for a magnetic field: 



/ 2{uj2roy - ujixz) \ 
-2{u2rox ^uiyz) 

\co;i(-rg +^2 - z‘^) J 


(16) 


This field is cylindrically symmetric as can be seen by 
the absence of a ^-dependence if equation is put in 
cylindrical coordinates: 


/ X 4rn\/a , ^ - 

^(-2w2ror(^ 

TT (r^ + r2 + 2:2) 

— 2LOizrr + LOi{—rQ + — z^)z). (17) 

Every field line lies on the surface of a member of a set 
of nested tori. The smallest reduces to a circle with radius 
of ro (magnetic axis), and the largest is a line along the 
^-axis. The field lines wind around the poloidal direction 
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with frequency uJi , and toroidal direction with frequency 
co’ 2 . If — is rational, ^ ^ and every field line is a 

(n, m) torus knot. We stress that every integral curve of 
this field is itself a knot (or circle, or ergodically spanning 
the toroidal surface), but the global field is smooth and 
continuous. The ratio ^ gives the ratio of toroidal to 
poloidal winding of the curves or the rotational transform 

i. 



CALCULATING THE HELICITY AND 
MAGNETIC ENERGY OF THE SAGDEEV FIELD 

The vector potential corresponding with the field in 
equation [Tbl defined as V x A = B, is given by 




A = 




2{rQL0iy - lo2Xz) 
-2{roU)ix + co2yz) 


7r(ro + r ) ^ 2,2 _|_ ^2 _ ^ 2 ^^ 

The inner product A • B is given by 


(18) 


A B = 


4:arluJiUJ2 

7r^(rQ + 


The helicity is then 


^ f 4arJ 

J 


The value of is given by 

16ar, 


4arlujiL02 ,3 4 

da; = arQU!iC02- 


(19) 


( 20 ) 


FIG. 8. fitting routine for finding the smallest torus. The 
magenta arrow indicates the vector fi, and the cyan arrow 
indicates the vector rop, which is the seed point for the next 
iteration in the fitting routine for the torus. 


An approximate major radius of the torus can be found 
by calculating 


ro = (xi - xo). 


(24) 


B B = 


222 ( 2-2 + b )6 ( 4 ^ 0(^^2 “ + V ^) + 

( 21 ) 


To find the orientation of the torus we calculate the cross 
product between the vector from the center to a point on 
the torus (x^ — xq) and a vector from a point on the torus 
to the next point (x^+i — x^) to get a vector pointing out 
of the torus and average this over all points on the field 
line: 


which allows us to calculate the integral of over all 
space 


/■ 


d'^x = aro(co’i + 


( 22 ) 


FINDING THE SMALLEST INVARIANT TORUS 
AND O RIE NTATION OF THE MAGNETIC 
STRUCTURE 


((Xj - Xq) X (Xj+i - Xj)) 
|(Xi - Xo) X (Xj+i - Xj)| ' 


(25) 


A vector perpendicular to the normal vector is found by 
solving 



(26) 


for c and defining the vector 


Properties of the field structure are extracted from the 
simulations. First the radius of the magnetic axis is found 
by a fitting routine. 

A point near the magnetic axis is guessed, from that 
point a the field is integrated using a fixed-step size field 
integration routine for a distance of 500/o- A set of points 
Xi on this field line is returned. The geometrical center 
of the torus is easily calculated by the weighted average 
of all the points: 


P = 


V2V^ 


(27) 


This yields a perpendicular vector in all cases except if 

h\ 

the normal vector is exactly in the 1 direction. 

Using the properties calculated above, a new starting 
point is taken to be 


Xo = (Xi). 


(23) 


s = Xo + rop. 


(28) 
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This starting point lies inside of the volume spanned by 
the previous torus. Integrating a new field line from this 
point yields a new field torus, but since the starting point 
is inside the previous torus this torus is nested inside the 
previous. This routine is repeated until the difference 
between two subsequent starting points reaches below a 
threshold of 0.005/o- 

Figure shows the first few steps of this iterative pro¬ 
cess, with n and s starting from xq drawn in magenta 
respectively in cyan. These are the properties extracted 
from the simulation: the radius of the smallest invari¬ 
ant torus given by |s — xo|, the geometrical center of the 
invariant torus xq and it’s orientation given by h 


FINDING THE ROTATIONAL TRANSFORM OF 
A SURFACE AND PARAMETERS FOR 
ANALYTICAL EXPRESSION 


of lOOO/o- The rotational transform is roughly constant 
on the magnetic surfaces, so it is sufficient to extract it 
only from a single curve. The rotational transform i is 
then equal to the ratio of poloidal to toroidal winding, 
and can be used to set the ratio ofco’i/co’ 2 . 


The value of ro is the radius of the smallest of the 
nested tori, and is found by using the fitting routine de¬ 
scribed above. The value of / d^x is calculated for 
the entire field in the simulation, and that value is used 
as the answer for equation 22 A field approximating the 
output of the simulation will have values of uji and UJ 2 
that are then uniquely defined through the ratio ^ and 
equation and the sign of the total helicity in equation 
3 of the paper. The results for all simulations are given 
in table [Til 


The simulation with n = 2 and T = 0 did not become 
ordered enough to be analyzed with the described rou- 


The rotational transform of field lines spanning invari¬ 
ant tori is found in the following way: First we calculate 
the number of times the field line crosses the plane de¬ 
fined by Xq and h. This number is twice the the toroidal 
winding rip, and is found by calculating the distance of 
each point in the curve to this plane and counting the 
zero crossings. The number of times the curve crosses the 
plane defined by xq and p is twice the toroidal winding 
nt. For a sufficiently long field line the ratio of poloidal 
to toroidal winding approaches the rotational transform 

rip/n, ~ i. 

We describe how to calculate the variables cji, UJ 2 and 
ro of equation [T6| of the paper, that describe a field with 
the same helicity, magnetic energy, and rotational trans¬ 
form of the field lines as the simulation field. In the 
dimensionless units of our simulations we can take the 
constant ^/a to be equal to I. 

To find the ratio of toroidal winding to poloidal wind¬ 
ing, we calculate the rotational transform of a single field 
line with a starting point St = xq + I.ISroP and a length 


TABLE II. Field Properties 


n 

T 

t 

ro 

uJi 

UJ2 

^ 0(^1 + ^ 2 ) 

UJ2 

2 

0 

- 

- 

- 

- 

- 

- 

3 

0 

54.0tA 

1.25 

0.049 

0.090 

0.020 

0.54 

4 

0 

22Ma 

1.10 

0.17 

0.19 

0.085 

0.90 

5 

0 

22Ma 

1.13 

0.22 

0.23 

0.I4I 

0.96 

6 

0 

22Ma 

1.14 

0.26 

0.25 

0.I9I 

1.05 

3 

0.9 

22Ma 

0.90 

0.18 

0.19 

0.048 

0.95 

3 

1.8 

22Ma 

0.78 

0.24 

0.27 

0.063 

0.89 

3 

2.6 

22Ma 

0.74 

0.30 

0.32 

0.075 

0.92 

3 

3.5 

22Ma 

0.75 

0.32 

0.33 

0.086 

0.98 

3 

4.4 

22Ma 

0.80 

0.32 

0.30 

0.096 

1.05 


tine before the simulation was stopped at t = GO^a, the 
simulation with n = 3 and T = 0 only at t = 54tA- This 
explains the different value for ro and much smaller value 
for ujiluj 2 ‘ All simulations show the emergence of a con¬ 
figuration where ro is roughly I, and ro becomes smaller 
if T is larger. The value of ujiIuj 2 is also around one, and 
increases with higher initial helicity. 















